#Function to print correlation between two covariates over time (GSS)
OverTimeCorr<-function(cov1, cov2){
  dd<-data%>%select_(cov1, cov2, "year")%>%filter(complete.cases(.))
  year<-names(table(dd$year))  #Summarize years from eligible data
  corr<-rep(NA, length(year))
  ##Couldn't get this to work with tapply
  for(i in 1:length(year)){
    corr[i]<-cor(dd[dd$year==year[i],1], dd[dd$year==year[i],2])
  }
  dd<-cbind(as.numeric(year), as.numeric(corr))
  dd<-as.data.frame(dd)
  dd[,1]<-as.numeric(as.character(dd[,1]))
  dd[,2]<-as.numeric(as.character(dd[,2]))
  dd<-na.omit(dd)
  names(dd)<-c("year", "corr")
  return(dd)
} 

### Function to code zero to 1 ####
zero.one<-function(x){
  min.x<-min(x, na.rm=T)
  max.x<-max(x-min.x, na.rm=T)
  return((x-min.x)/max.x)
}
##### ANES data. Final #####
# Marginal mlogit results, For Figure 1
sm.anes.margins<-function(output, sm.value=0){
  t<-model.matrix(output)  # The data
  o<-cbind(
    1,
    sm.value, #SM
    1,
    sm.value*1, # SM x Conservative
    t[,5:dim(t)[2]]
  )
  p<-cbind(
    1,
    sm.value, #SM
    0,
    sm.value*0, # SM x Conservative
    t[,5:dim(t)[2]]
  )
    beta.sim<-mvrnorm(1000, c(coef(output)[1,], coef(output)[2,]), vcov(output)) ##Draw samples from multivariate distrbution
    ###Generate the Predictions###
    g1.dem.o<-o%*%t(beta.sim[,1:11])
    g1.dem.p<-p%*%t(beta.sim[,1:11])
    g1.ind.o<-o%*%t(beta.sim[,12:22])
    g1.ind.p<-p%*%t(beta.sim[,12:22])
# Generate averaged estimates     
    pr.dem.o<-apply(exp(g1.dem.o)/(1+exp(g1.ind.o)+exp(g1.dem.o)), 2, mean)
    pr.dem.p<-apply(exp(g1.dem.p)/(1+exp(g1.ind.p)+exp(g1.dem.p)), 2, mean)
    pr.ind.o<-apply(exp(g1.ind.o)/(1+exp(g1.ind.o)+exp(g1.dem.o)), 2, mean)
    pr.ind.p<-apply(exp(g1.ind.p)/(1+exp(g1.ind.p)+exp(g1.dem.p)), 2, mean)
    pr.rep.o<-apply((1)/(1+exp(g1.ind.o)+exp(g1.dem.o)), 2, mean)
    pr.rep.p<-apply((1)/(1+exp(g1.ind.p)+exp(g1.dem.p)), 2, mean)
## Generate marginal effects     
    g1<-pr.dem.o-pr.dem.p
    g2<-pr.ind.o-pr.ind.p
    g3<-pr.rep.o-pr.rep.p
## Generate list of responses from this
# Key: 
        return(list(dem.margin=data.frame(min.2.5=quantile(g1, 0.025),
                                          min.25=quantile(g1, 0.25),
                                          max.75=quantile(g1, 0.75),
                                          max.97.5=quantile(g1, 0.975),
                                          mean.score=quantile(g1, 0.5)),
                    ind.margin=data.frame(min.2.5=quantile(g2, 0.025),
                                          min.25=quantile(g2, 0.25),
                                          max.75=quantile(g2, 0.75),
                                          max.97.5=quantile(g2, 0.975),
                                          mean.score=quantile(g2, 0.5)),
                    rep.margin=data.frame(min.2.5=quantile(g3, 0.025),
                                          min.25=quantile(g3, 0.25),
                                          max.75=quantile(g3, 0.75),
                                          max.97.5=quantile(g3, 0.975),
                                          mean.score=quantile(g3, 0.5)),
                    prediction.demhigh=data.frame(min.2.5=quantile(pr.dem.o, 0.025),
                                          min.25=quantile(pr.dem.o, 0.25),
                                          max.75=quantile(pr.dem.o, 0.75),
                                          max.97.5=quantile(pr.dem.o, 0.975),
                                          mean.score=quantile(pr.dem.o, 0.5)),
                    prediction.demlow=data.frame(min.2.5=quantile(pr.dem.p, 0.025),
                                                  min.25=quantile(pr.dem.p, 0.25),
                                                  max.75=quantile(pr.dem.p, 0.75),
                                                  max.97.5=quantile(pr.dem.p, 0.975),
                                                  mean.score=quantile(pr.dem.p, 0.5)),
                    prediction.rephigh=data.frame(min.2.5=quantile(pr.rep.o, 0.025),
                                                  min.25=quantile(pr.rep.o, 0.25),
                                                  max.75=quantile(pr.rep.o, 0.75),
                                                  max.97.5=quantile(pr.rep.o, 0.975),
                                                  mean.score=quantile(pr.rep.o, 0.5)),
                    prediction.replow=data.frame(min.2.5=quantile(pr.rep.p, 0.025),
                                                  min.25=quantile(pr.rep.p, 0.25),
                                                  max.75=quantile(pr.rep.p, 0.75),
                                                  max.97.5=quantile(pr.rep.p, 0.975),
                                                  mean.score=quantile(pr.rep.p, 0.5))
 ))
    
  }
### Voting ####
sm.anes.margins2<-function(output, sm.value=0){
  t<-model.matrix(output)  # The data
  o<-cbind(
    1,
    sm.value, #SM
    1,
    sm.value*1, # SM x Conservative
    t[,5:dim(t)[2]]
  )
  p<-cbind(
    1,
    sm.value, #SM
    0,
    sm.value*0, # SM x Conservative
    t[,5:dim(t)[2]]
  )
  beta.sim<-mvrnorm(1000, coef(output), vcov(output)) ##Draw samples from multivariate distrbution
  ###Generate the Predictions###
  pr.rep.o<-apply(plogis(o%*%t(beta.sim)), 2, mean)
  pr.rep.p<-apply(plogis(p%*%t(beta.sim)), 2, mean)
  ## Generate marginal effects     
  g1<-pr.rep.o-pr.rep.p
  ## Generate list of responses from this
  # Key: 
  return(list(margins=data.frame(min.2.5=quantile(g1, 0.025),
                                    min.25=quantile(g1, 0.25),
                                    max.75=quantile(g1, 0.75),
                                    max.97.5=quantile(g1, 0.975),
                                    mean.score=quantile(g1, 0.5)),
              prediction.rephigh=data.frame(min.2.5=quantile(pr.rep.o, 0.025),
                                            min.25=quantile(pr.rep.o, 0.25),
                                            max.75=quantile(pr.rep.o, 0.75),
                                            max.97.5=quantile(pr.rep.o, 0.975),
                                            mean.score=quantile(pr.rep.o, 0.5)),
              prediction.replow=data.frame(min.2.5=quantile(pr.rep.p, 0.025),
                                           min.25=quantile(pr.rep.p, 0.25),
                                           max.75=quantile(pr.rep.p, 0.75),
                                           max.97.5=quantile(pr.rep.p, 0.975),
                                           mean.score=quantile(pr.rep.p, 0.5))
  ))
  
}
### polarization ####
sm.anes.margins3<-function(output, sm.value=0){
  t<-model.matrix(output)  # The data
  o<-cbind(
    1,
    sm.value, #SM
    1,
    sm.value*1, # SM x Conservative
    t[,5:dim(t)[2]]
  )
  p<-cbind(
    1,
    sm.value, #SM
    min(t[,2]),
    sm.value*0, # SM x Conservative
    t[,5:dim(t)[2]]
  )
  beta.sim<-mvrnorm(1000, coef(output), vcov(output)) ##Draw samples from multivariate distrbution
  ###Generate the Predictions###
  pr.rep.o<-apply((o%*%t(beta.sim)), 2, mean)
  pr.rep.p<-apply((p%*%t(beta.sim)), 2, mean)
  ## Generate marginal effects     
  g1<-pr.rep.o-pr.rep.p
  ## Generate list of responses from this
  # Key: 
  return(list(margins=data.frame(min.2.5=quantile(g1, 0.025),
                                 min.25=quantile(g1, 0.25),
                                 max.75=quantile(g1, 0.75),
                                 max.97.5=quantile(g1, 0.975),
                                 mean.score=quantile(g1, 0.5)),
              prediction.rephigh=data.frame(min.2.5=quantile(pr.rep.o, 0.025),
                                            min.25=quantile(pr.rep.o, 0.25),
                                            max.75=quantile(pr.rep.o, 0.75),
                                            max.97.5=quantile(pr.rep.o, 0.975),
                                            mean.score=quantile(pr.rep.o, 0.5)),
              prediction.replow=data.frame(min.2.5=quantile(pr.rep.p, 0.025),
                                           min.25=quantile(pr.rep.p, 0.25),
                                           max.75=quantile(pr.rep.p, 0.75),
                                           max.97.5=quantile(pr.rep.p, 0.975),
                                           mean.score=quantile(pr.rep.p, 0.5))
  ))
  
}
### Samara Klar's 2014 Data #####
treatment.effect<-function(output,  sm){
  require(MASS)
  t<-model.matrix(output)[,-1]
  homo<-cbind(
    sm,
    0,
    1,
    sm*0,
    sm*1
  )
  hetero<-cbind(
    sm,
    1,
    0,
    sm*1,
    sm*0
  )
  control<-cbind(
    sm,
    0,
    0,
    sm*0,
    sm*0
  )
  set.seed<-123
  beta.sim<-mvrnorm(1000, c(output$coefficients, output$zeta), vcov(output)) ##Draw samples from multivariate distrbution
  g1a<-plogis(cbind(homo, -1, 0, 0, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g2a<-plogis(cbind(homo, 0, -1, 0, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g3a<-plogis(cbind(homo, 0, 0, -1, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g4a<-plogis(cbind(homo, 0, 0, 0,-1, 0)%*%t(beta.sim))  ### The average effect. Simulate b and auth
  g5a<-plogis(cbind(homo, 0, 0, 0, 0, -1)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  
  
  g1b<-plogis(cbind(hetero, -1, 0, 0, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g2b<-plogis(cbind(hetero, 0, -1, 0, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g3b<-plogis(cbind(hetero, 0, 0, -1, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g4b<-plogis(cbind(hetero, 0, 0, 0,-1, 0)%*%t(beta.sim))  ### The average effect. Simulate b and auth
  g5b<-plogis(cbind(hetero, 0, 0, 0, 0, -1)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  
  g1c<-plogis(cbind(control, -1, 0, 0, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g2c<-plogis(cbind(control, 0, -1, 0, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g3c<-plogis(cbind(control, 0, 0, -1, 0, 0)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  g4c<-plogis(cbind(control, 0, 0, 0,-1, 0) %*%t(beta.sim))  ### The average effect. Simulate b and auth
  g5c<-plogis(cbind(control, 0, 0, 0, 0, -1)%*%t(beta.sim)) ### The average effect. Simulate b and auth
  
  pr.1<-g4a
  pr.2<-g4b  ###  
  pr.3<-g4c  ### Prioritize Oil energy 
  
  homo.treat<-pr.1-pr.3
  hetero.treat<-pr.2-pr.3
  
  
  return(list(homo=data.frame(min.2.5=quantile(pr.1, 0.025),
                              min.25=quantile(pr.1, 0.25),
                              max.75= quantile(pr.1, 0.75),
                              max.97.5=quantile(pr.1, 0.975),
                              mean.score=quantile(pr.1, 0.5)),
              hetero=data.frame(min.2.5=quantile(pr.2, 0.025),
                                min.25=quantile(pr.2, 0.25),
                                max.75= quantile(pr.2, 0.75),
                                max.97.5=quantile(pr.2, 0.975),
                                mean.score=quantile(pr.2, 0.5)),
              control=data.frame(min.2.5=quantile(pr.3, 0.025),
                                 min.25=quantile(pr.3, 0.25),
                                 max.75= quantile(pr.3, 0.75),
                                 max.97.5=quantile(pr.3, 0.975),
                                 mean.score=quantile(pr.3, 0.5)),
              homo.treat=data.frame(min.2.5=quantile(homo.treat, 0.025),
                                    min.25=quantile(homo.treat, 0.25),
                                    max.75= quantile(homo.treat, 0.75),
                                    max.97.5=quantile(homo.treat, 0.975),
                                    mean.score=quantile(homo.treat, 0.5)),
              hetero.treat=data.frame(min.2.5=quantile(hetero.treat, 0.025),
                                      min.25=quantile(hetero.treat, 0.25),
                                      max.75= quantile(hetero.treat, 0.75),
                                      max.97.5=quantile(hetero.treat, 0.975),
                                      mean.score=quantile(hetero.treat, 0.5))
  ))
}
##CCAP Functions###
pred.generator<-function(output, model="ologit", democrat=1, independent=1){
  require(MASS)
  if(model=="ologit"){
    t<-model.matrix(output)[,-1]
    p<-cbind(
      1,
      democrat,
      independent,
      democrat*1,
      independent*1,
      t[,6:dim(t)[2]])
    o<-cbind(
      0,
      democrat,
      independent,
      democrat*0,
      independent*0,
      t[,6:dim(t)[2]])
    set.seed<-123
    beta.sim<-mvrnorm(1000, c(output$coefficients, output$zeta), vcov(output)) ##Draw samples from multivariate distrbution
    g1<-apply(plogis(cbind(p, -1, 0, 0, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g2<-apply(plogis(cbind(p, 0, -1, 0, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g3<-apply(plogis(cbind(p, 0, 0, -1, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g4<-apply(plogis(cbind(p, 0, 0, 0,-1, 0, 0)%*%t(beta.sim)), 2, mean)  ### The average effect. Simulate b and auth
    g5<-apply(plogis(cbind(p, 0, 0, 0, 0, -1, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g6<-apply(plogis(cbind(p, 0, 0, 0, 0, 0, -1)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g1a<-apply(plogis(cbind(o, -1, 0, 0, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g2a<-apply(plogis(cbind(o, 0, -1, 0, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g3a<-apply(plogis(cbind(o, 0, 0, -1, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g4a<-apply(plogis(cbind(o, 0, 0, 0,-1, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g5a<-apply(plogis(cbind(o, 0, 0, 0, 0, -1, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g6a<-apply(plogis(cbind(o, 0, 0, 0, 0, 0, -1)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    pr.1<-1-g2
    pr.2<-1-g2a  ### Happy 
    marg1<-pr.1-pr.2
    pr.3<-g5
    pr.4<-g5a ### Unhappy
    marg2<-pr.3-pr.4
    
    return(list(Happy.High=data.frame(min.2.5=quantile(pr.1, 0.025),
                                      min.25=quantile(pr.1, 0.25),
                                      max.75= quantile(pr.1, 0.75),
                                      max.97.5=quantile(pr.1, 0.975),
                                      mean.score=quantile(pr.1, 0.5)),
                Happy.Low=data.frame(min.2.5=quantile(pr.2, 0.025),
                                     min.25=quantile(pr.2, 0.25),
                                     max.75= quantile(pr.2, 0.75),
                                     max.97.5=quantile(pr.2, 0.975),
                                     mean.score=quantile(pr.2, 0.5)),
                Unhappy.High=data.frame(min.2.5=quantile(pr.3, 0.025),
                                        min.25=quantile(pr.3, 0.25),
                                        max.75= quantile(pr.3, 0.75),
                                        max.97.5=quantile(pr.3, 0.975),
                                        mean.score=quantile(pr.3, 0.5)),
                Unhappy.Low=data.frame(min.2.5=quantile(pr.4, 0.025),
                                       min.25=quantile(pr.4, 0.25),
                                       max.75= quantile(pr.4, 0.75),
                                       max.97.5=quantile(pr.4, 0.975),
                                       mean.score=quantile(pr.4, 0.5)),
                marginal.happy=data.frame(min.2.5=quantile(marg1, 0.025),
                                          min.25=quantile(marg1, 0.25),
                                          max.75= quantile(marg1, 0.75),
                                          max.97.5=quantile(marg1, 0.975),
                                          mean.score=quantile(marg1, 0.5)),
                marginal.sad=data.frame(min.2.5=quantile(marg2, 0.025),
                                        min.25=quantile(marg2, 0.25),
                                        max.75= quantile(marg2, 0.75),
                                        max.97.5=quantile(marg2, 0.975),
                                        mean.score=quantile(marg2, 0.5))
    ))
  }
  if(model=="logit"){
    t<-model.matrix(output)
    p<-cbind(
      1,
      1,
      democrat,
      independent,
      democrat*1,
      independent*1,
      t[,7:dim(t)[2]])
    o<-cbind(
      1,
      0,
      democrat,
      independent,
      democrat*0,
      independent*0,
      t[,7:dim(t)[2]])
    set.seed<-123
    beta.sim<-mvrnorm(1000, output$coefficients, vcov(output)) ##Draw samples from multivariate distrbution
    g1<-apply(plogis(p%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g2<-apply(plogis(o%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    me<-g1-g2
    return(list(High=data.frame(min.2.5=quantile(g1, 0.025),
                                min.25=quantile(g1, 0.25),
                                max.75= quantile(g1, 0.75),
                                max.97.5=quantile(g1, 0.975),
                                mean.score=quantile(g1, 0.5)),
                Low=data.frame(min.2.5=quantile(g2, 0.025),
                               min.25=quantile(g2, 0.25),
                               max.75= quantile(g2, 0.75),
                               max.97.5=quantile(g2, 0.975),
                               mean.score=quantile(g2, 0.5)),
                marginal=data.frame(min.2.5=quantile(me, 0.025),
                                    min.25=quantile(me, 0.25),
                                    max.75= quantile(me, 0.75),
                                    max.97.5=quantile(me, 0.975),
                                    mean.score=quantile(me, 0.5))
                
    ))
  }
}
#### Reported Neighborhood ####
pred.generator.friends<-function(output, model="ologit",  democrat=1, independent=1){
  require(MASS)
  if(model=="ologit"){
    t<-model.matrix(output)[,-1]
    p<-cbind(
      1,
      democrat,
      independent,
      democrat*1,
      independent*1,
      t[,6:dim(t)[2]])
    o<-cbind(
      0,
      democrat,
      independent,
      democrat*0,
      independent*0,
      t[,6:dim(t)[2]])
    set.seed<-123
    beta.sim<-mvrnorm(1000, c(output$coefficients, output$zeta), vcov(output)) ##Draw samples from multivariate distrbution
    g1<-apply(plogis(cbind(p, -1, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g2<-apply(plogis(cbind(p, 0, -1, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g3<-apply(plogis(cbind(p, 0, 0, -1, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g4<-apply(plogis(cbind(p, 0, 0, 0,-1)%*%t(beta.sim)), 2, mean)  ### The average effect. Simulate b and auth
    
    g1a<-apply(plogis(cbind(o, -1, 0, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g2a<-apply(plogis(cbind(o, 0, -1, 0, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g3a<-apply(plogis(cbind(o, 0, 0, -1, 0)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    g4a<-apply(plogis(cbind(o, 0, 0, 0,-1)%*%t(beta.sim)), 2, mean) ### The average effect. Simulate b and auth
    
    pr.1<-1-g2
    pr.2<-1-g2a  ### Happy 
    marg1<-pr.1-pr.2
    pr.3<-g3
    pr.4<-g3a ### Unhappy
    marg2<-pr.3-pr.4
    return(list(Left.High=data.frame(min.2.5=quantile(pr.1, 0.025),
                                     min.25=quantile(pr.1, 0.25),
                                     max.75= quantile(pr.1, 0.75),
                                     max.97.5=quantile(pr.1, 0.975),
                                     mean.score=quantile(pr.1, 0.5)),
                Left.Low=data.frame(min.2.5=quantile(pr.2, 0.025),
                                    min.25=quantile(pr.2, 0.25),
                                    max.75= quantile(pr.2, 0.75),
                                    max.97.5=quantile(pr.2, 0.975),
                                    mean.score=quantile(pr.2, 0.5)),
                Right.High=data.frame(min.2.5=quantile(pr.3, 0.025),
                                      min.25=quantile(pr.3, 0.25),
                                      max.75= quantile(pr.3, 0.75),
                                      max.97.5=quantile(pr.3, 0.975),
                                      mean.score=quantile(pr.3, 0.5)),
                Right.Low=data.frame(min.2.5=quantile(pr.4, 0.025),
                                     min.25=quantile(pr.4, 0.25),
                                     max.75= quantile(pr.4, 0.75),
                                     max.97.5=quantile(pr.4, 0.975),
                                     mean.score=quantile(pr.4, 0.5)),
                marginal.left=data.frame(min.2.5=quantile(marg1, 0.025),
                                         min.25=quantile(marg1, 0.25),
                                         max.75= quantile(marg1, 0.75),
                                         max.97.5=quantile(marg1, 0.975),
                                         mean.score=quantile(marg1, 0.5)),
                marginal.right=data.frame(min.2.5=quantile(marg2, 0.025),
                                          min.25=quantile(marg2, 0.25),
                                          max.75= quantile(marg2, 0.75),
                                          max.97.5=quantile(marg2, 0.975),
                                          mean.score=quantile(marg2, 0.5))
                
                
                
    ))
  }
}
